This is a working protocol for processing terrestrial lidar scans that we collect. By working protocol I mean that all aspects of the process can be improved! We can all find ways to amend and improve it. Please speak up if you have suggestions on how to make the protocol more user-friendly, more efficient, or more accurate! And also let us know if you run into issues.
Overview
Step-by-step
Is there a way to judge particularly bad scans Tanner?
The terrestrial lidar scans we collect have a few issues we need to deal with.
First, scans can include lidar point returns up to 400 m, but only the first ~60 m of returns is usable. After that, returns are too sparse. Further, small errors in alignment (like a slight tilt) are magnified at larger distances
Raw lidar scan showing rapid decay of return density
Second, lidar scans are taken from a fixed point and a single sscan can only see one side of an object. Multiple scans must be combined to capture all sides of objects.
View beneath a tree from a single scan showing C shaped pattern of returns from one side of the tree
Third, lidar scans are huge. A small area (0.16 ha or 40 \(\times\) 40 m area) typically comprises nearly 1GB of data. Proccessing this takes a lot of computer memory, so we must break the data into tiles (or lidR::chunks) to process a bit at a time.
Forth, when dealing with lidar data there are often artefacts to deal with including noise points introduced during collection, and misregistration errors that can occur.
Each of these issues can be dealt with by pre-processing the scans with the following steps.
We’ll mostly be using R Statistical software for processing lidar data. The functions in R come from small bits of code (or packages) contributed by other R users. You will need to download and install these packages once on each computer you use for processing data. We will be downloading and installing some “unauthorized” R packages, and for that you’ll need to have Rtools installed on your computer
RtoolsRtoolsPATH="${RTOOLS40_HOME}\usr\bin;${PATH}".Renviron in your Documents folder.R packagesOnce installed you can install the R packages you’ll need using the following code:
# Run this once on each computer
# main package for dealing with lidar data in R
install.packages('lidR', repos="https://cloud.r-project.org")
# fast software for handling GIS data
install.packages('sf', repos="https://cloud.r-project.org")
# software for installing `unauthorized` R pakcages
install.packages('devtools', repos="https://cloud.r-project.org")
The remaining packages are not on CRAN, but hosted by users on github pages. The devtools package allows us to install these directly.
# Forestry specific Lidar functions for TLS
devtools::install_github('tiagodc/TreeLS')
# good tree detection software (our modified version of spanner)
devtools::install_github('jbcannon/spanner')
# functions for our labs common tasks
devtools::install_github('jbcannon/landecoutils')
To merge, clip, and tile scans from a large plot (> 0.2 ha) use the following R Code
ctg Path to a folder containing all *.las files for a plot to be combinedout_dirbnd A shapefile of the plot boundary. Note in this example we subset to a specific plottile_size side length of tiles that the data is cut intobuffer width of a buffer around the plot boundary bnd. This is included to reduce edge effects (e.g., overhanging crowns). By default buffer=10 but is an optional argumentmax_scan_distance maximum radius from scanner within which lidar returns are included. Defaults to max_scan_distance = 60 which is working for us.scan_locations this is an optional input of an sf object containing the coordinates for all scans. This is useful for quickly identifying for a particular region, which scans are within the max_scan_distance for inclusion. Scan locations can be passed from landecoutils::find_ctg_centroids. If NULL the scan locations will be computed within the stitch_TLS... function. It’s just provided separately to increase effeciency if the same ctg is used more than once.# load required packages
library(lidR)
library(sf)
library(landecoutils)
#bring up the help file for the main function used...
?stitch_TLS_dir_to_LAS_tiles
?find_ctg_centroids
# Path to folder containing the scans
ctg = readLAScatalog('E:/ecofor_segmentation/raw_scans/Plot_3')
# shapefile indicating Plot boundary. Here we have to subet Plot 3
bnd = st_read(
'R:/landscape_ecology/projects/eco_for/overstory_pattern/shp/ecoforplots.shp')
bnd = subset(bnd, PlotID==3)
# This function finds the "center" of all scans so we can delete points > 60 m away
scan_locations = find_ctg_centroids(ctg)
stitch_TLS_dir_to_LAS_tiles(ctg = ctg,
out_dir = 'E:/ecofor_segmentation/tiled_scans/Plot_3',
bnd = bnd,
tile_size = 30,
buffer = 10, #defaults to 10
max_scan_distance = 60, # defaults to 60
scan_locations = scan_locations)
Demonstration of stitch_TLS_dir_to_LAS_tiles function showing scan locations (blue), over study boundary (thick rectange), and tiled grid of study area.
The function works by automatically setting up a tiled grid over the study boundary (bnd) and buffer you provide as specified by tile_size. For each tile in the grid, the function loads and clips any scans overlapping the tile (blue circles) and merges them together. Depending on the scanner used, some tiles have different column names and the function rectifies that by choosing data common to all scans. Each tile is saved in the directory specified by out_dir with file name basedon on the XY coordinate of the bottom left corner of each tile.
Example output of stitch_TLS_dir_to_LAS_tiles with 30 m tile size
For smaller areas (e.g., a fixed 18 m radius plot) you may be able to avoid tiling the area, and just merge together multiple scans into a single *.las file. For this use the similar landecoutils::stitch_TLS_dir_to_LAS.
Inputs are the same as above for stitch_TLS_dir_to_LAS_tiles except that the plot boundary or region of interest is defined by roi rather than bnd and there is no need to specify a tile_size. A small buffer (usually 10 m) is still specified to account for edge effects.
library(lidR)
library(sf)
library(landecoutils)
#bring up the help file for the main functions used...
?stitch_TLS_dir_to_LAS
?find_ctg_centroids
# Path to folder containing the scans of a small plot area
ctg = readLAScatalog('E:/ecofor_segmentation/raw_scans/Plot_11') # Path to folder containing the scans
# Load a shapefile with the area of interest (e.g., ~0.2 ha)
bnd = st_read('R:/landscape_ecology/projects/eco_for/click_mapping_2021/click_mapping_ground_truth/grid-cells.shp')
bnd = subset(bnd, Plot==11 & cell == 90) #subsets bnd to only the Plot and cell of interest.
scan_locations = find_ctg_centroids(ctg)
stitch_TLS_dir_to_LAS(ctg = ctg,
out_las = 'E:/ecofor_segmentation/ground_truth_scans/Plot11_cell90.las', #where to write the .las
roi = bnd,
buffer = 10,
max_scan_distance = 60,
scan_locations = scan_locations)
Some scans will include meaningless returns (noise) that are easily noticeable when viewed in CloudCompare with an increased point size. These noisy returns come from smoke, water vapor, or scan detector issues and will need to be removed. Lidar returns with very low intensity (e.g., \(< 4000\)) are often noise and can be removed. We also use a noise filter found in the lidR::classify_noise() function. It uses an isolated voxel filter (ivf()) to remove all points that don???t have a certain number of points (n) in the neighboring voxels with size specified by res. I have found that ivf(res = 0.1, n = 15) works for most issues encountered to date.
Example of tiled TLS scan with noise to be removed
library(lidR)
# Code to do noise removal after merging (especially 2020 TallTimbers scans)
las_fn = 'E:/ecofor_segmentation/ground_truth_scans/Plot11_cell90.las'
new_fn = 'E:/ecofor_segmentation/ground_truth_scans/Plot11_cell90_denoised.las'
las = readLAS(las_fn)
las = filter_poi(las, Intensity < 4000) #18 is the LAS specification for noise
las = classify_noise(las, ivf(res=0.1,n=15)) #This setting seems to work for our TT scans with noise
las = filter_poi(las, Classification !=18) #18 is the LAS specification for noise
writeLAS(las, new_fn) #overwrite the old file with the filtered one.
You can remove noise from a single file at a time, or you can use a for loop in R to denoise an entire group at a time. Here is an example of de-noising all files in a folder with ???Plot6_??? in their filename. Replace pattern to pattern match your files of interest. Test out the ivf settings on a single file first as this code overwrites original tiles after denoising.
library(lidR)
fn = list.files(out_dir, pattern='Plot6_', full.names=TRUE)
fn = grep('.las', fn, value=TRUE)
for(f in fn) {
cat('working on ', basename(f),'...\n')
x = lidR::filter_poi(lidR::readLAS(f), Intensity>4000)
x = lidR::filter_poi(lidR::classify_noise(x, lidR::ivf(res=0.1, n=15)),
Classification!=18)
lidR::writeLAS(x, f, index=TRUE)
}
Because each area of interest is composed of two or more scans from different angles, sometimes not everything lines up perfectly. If a scan is misaligned it can lead to trees looking much thicker than they really are, or if more severe, lead to the duplication of some trees! Dealing with misregistration generally means removing the scan that is misaligned from the other merged scans. This may be easier done at the earlier stage of stitching scans in RiScan. However, they can also be removed later using the following procedure and code
Example of misregistration errors leading to duplicate trees
Note that not all mis-alignments are due to registration. Wind can move plant parts between scans and make it appear as though the scan has minor misregistation. In the following screenshot you can see that small branches appear to be misaligned, but large tree bowls do not. Another clue that this is from wind is that the ???phantom??? branches come from two different scans (blue and gold) so it is likely the branch moved between the scans.
The effect of this kind of misregistration is minimal, as we are generally focused on larger trees and trunks which don’t move in the wind. So you can ignore this type of registration error.
Misregistration due to wind movement. Note that trunk (white error) has no issues but small branches (red arrows) are apparently misaligned
The basic idea here is to find misregistered portions of scans, then isolate and remove the points coming from the offending scan. Here we identify bad portions of scans in CloudCompare, then remove them with R (though you can also remove in CC).
In CloudCompare find the areas of the scan that have misalginment. Set the Scalar Field to gpstime which logs the time at which each return was captured. Generally returns from a single scan are all grouped in within a short time. Adjust the display ranges of the gpstime attribute until you’ve identified the returns you want to remove. In the case below the green points need to be removed. These are prepresented by 3 peaks in the display range window.
Using the gpstime attribute to isolate duplicated and misregistered returnsin CloudCompare
Histogram of gpstime attribute and area to be filtered out
Now using R, you can load the tile, view a histogram of the gpstime attribute, note the range of times to remove, and filter them out using the following code as an example. Use the filter_poi function and create a gpstime function to keep only the appropriate range. Note the the values may not match numerically because CloudCompare rescales some attributes. So look at the histogram and find the range that shape that corresponds. Usually the time ranges are in short spikes representing when the scanner is running and breaks during transportation.
Within the filter_poi() function, you’ll need to write a logical statement to remove the correct pixels. & is the AND operator and | is the OR operator. See other example logical operations below
library(lidR)
x = readLAS('E:/ecofor_segmentation/ground_truth_scans/Plot5_cell7.las')
hist(x$gpstime, breaks=1e3)
x = filter_poi(x, gpstime < 400 | gpstime > 700)
writeLAS(x, 'E:/ecofor_segmentation/ground_truth_scans/Plot5_cell7.las', index=TRUE)
Examples of logical operators you may use.
gpstime below 700: gpstime < 700gpstime above 200: gpstime > 200gpstime outside of the range of 400-700: gpstime < 400 | gpstime > 700gpstime is outside fo 400-700 or outside of 100-200: (gpstime < 100 | gpstime > 200) & (gpstime < 400 | gpstime > 700)You have made it! Once scans are planned, collected, and pre-processed you have clean Lidar tiles that are ready to be processed. The stem mapping portion of the protocol should be considered the most tentative. It will more consistently change as we try new algorithms, and individual project goals change. But it is always useful to have a basic stem map for many of our scans.
Here is a summary of the main steps used for stem mapping, though it’s best to view the code directly to see the steps.
Diagram of cloth simulation filter (Figure from Zhang et al. 2016 Remote Sensing)
Example of intensity filter used to remove most vegetation and reveal tree boles
spanner does just that, but it also uses information on vertical alignment of points to help find trees. No algorithm is perfect. The settings below have done the best at capturing smaller trees.example of horizontal slice taken from 1 to 2 m to use for cylinder fitting in spanner
spanner package does this using graph theory. Trees are “grown” from their location on the ground and nearby points are successively added until all points are assigned.Example TLS scan with high sapling density segmented with spanner
Take look at the code below and see if you can identify where each of these steps take place. And find a few that I glossed over.
To stem map TLS tiles, use the following code. This works on small tiles, but you’ll need to use something like catalog apply for large areas.
install_github('tiagodc/TreeLS')
install_github('jbcannon/spanner') #slightly modified version of spanner
install_github('jbcannon/landecoutils')
library(lidR)
library(sf)
library(TreeLS)
library(landecoutils)
library(spanner)
?segment_with_spanner
bnd = sf::st_read('R:/landscape_ecology/projects/eco_for/click_mapping_2021/click_mapping_ground_truth/grid-cells.shp')
bnd = subset(bnd, Plot == 1 & cell == 46)
nCores=4
in_las = 'E:/ecofor_segmentation/ground_truth_scans/Plot1_cell46.las'
out_map = 'C:/Users/jeffery.cannon/OneDrive - Joseph W. Jones Ecological Research Center/Desktop/test.shp'
tree_seg = 'C:/Users/jeffery.cannon/OneDrive - Joseph W. Jones Ecological Research Center/Desktop/testseg.las'
landecoutils::segment_with_spanner(las = in_las, bnd=bnd, stemmap_shp = out_map, tree_seg_las = tree_seg)